//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

// try to see errors in bicubic interpolation for the cube polywell

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>

#include "ephi.hpp"

#define SPACING 0.01
#define RADIUS 0.15
#define WR 0.035

static prec_t maxbf_error, maxef_error;
static vect3d bferrpos, eferrpos;
bool do_relative;
bool discard_middle;

static void report_error(const char *name = "")
{
  char b1[128], b2[128];
  if ( do_relative )
    {
      printf ("%s: b-field %g%% %s, e-field %g%% %s\n",
          name, prec2double (maxbf_error), bferrpos.sprint(b1), prec2double (maxef_error), eferrpos.sprint(b2));
    }
  else
    {
      printf ("%s: b-field %g %s, e-field %g %s\n",
          name, prec2double (maxbf_error), bferrpos.sprint(b1), prec2double (maxef_error), eferrpos.sprint(b2));
    }
  maxbf_error = maxef_error = 0;
}

struct DiffStaticElement : public StaticElement
{
  DiffStaticElement(Statics& s1, Statics& s2) : s1(s1), s2(s2)
  {
  }

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield)
  {
    vect3d bf1, ef1, bf2, ef2;
    s1.getFields (pos, bf1, ef1);
    s2.getFields (pos, bf2, ef2);
    if ( do_relative )
      {
        if ( discard_middle && pos.length() < 0.025 ) // relative error in the middle can be big but its effect very small
            return;
        if ( bf1.length() != 0 )
            bfield += vect3d(0, 0, 1) * (bf1 - bf2).length() / bf1.length();
        if ( ef1.length() != 0 )
            efield += vect3d(0, 0, 1) * (ef1 - ef2).length() / ef1.length();

        if ( !inContact (pos) )
          {
            prec_t bferr = (bf1 - bf2).length() / bf1.length() * 100;
            if ( bferr > maxbf_error )
              {
                maxbf_error = bferr;
                bferrpos = pos;
              }
            prec_t eferr = (ef1 - ef2).length() / ef1.length() * 100;
            if ( eferr > maxef_error )
              {
                maxef_error = eferr;
                eferrpos = pos;
              }
          }
      }
    else
      {
        bfield += bf1 - bf2;
        efield += ef1 - ef2;

        if ( !inContact (pos) )
          {
            prec_t bferr = (bf1 - bf2).length();
            if ( bferr > maxbf_error )
              {
                maxbf_error = bferr;
                bferrpos = pos;
              }
            prec_t eferr = (ef1 - ef2).length();
            if ( eferr > maxef_error )
              {
                maxef_error = eferr;
                eferrpos = pos;
              }
          }
      }
  }

  void addPotential (const vect3d& pos, prec_t& potential)
  {
    potential += s1.getPotential(pos) - s2.getPotential(pos);
  }

  void draw(Screen& screen, const Screen::color& color)
  {
    s1.draw (screen, color);
  }

  std::string get_povray_script ()
  {
    return std::string();
  }

  bool inContact (const vect3d& pos) 
  {
    return s1.inContact (pos);
  }

  Statics &s1, &s2;
};

int main (int argc, char *argv[])
{
  Ephi::setDebugLevel (Ephi::DEBUG_APP);
  {
  CoilGen cgen;
  cgen.radius = RADIUS;
  cgen.wr = WR;
  cgen.do_efield = true;
  cgen.shape = CoilGen::ROUND;
  cgen.geom_radial = 4 * RADIUS;
  cgen.geom_axial = 4 * RADIUS;
  cgen.memory_usage = 500;
  cgen.coil_matrix_n = 1;
  cgen.do_cubic = true;
  cgen.generate ("test19_coil.cgen", false, false);

  CoilData cd;
  assert (cd.load ("test19_coil.cgen"));

  Statics s1;
  Statics s2;

  ShapedCoilFactory cf1;
  cf1.make (s1, vect3d(0, 0, 0), vect3d(1, 0, 0), RADIUS, WR, 1, 1);
  CoilDataFactory<StaticCICU> cf2(cd);
  cf2.make (s2, vect3d(0, 0, 0), vect3d(1, 0, 0), RADIUS, WR, 1, 1);

  Statics statics;
  ElectroDynamics dyn(statics);
  statics.addStaticElement (new DiffStaticElement (s1, s2));
  Scene scene(dyn);
  scene.set_coloring (0, true);

  Screen screen(800, 800, 0.7, 0.7);
  scene.calc (screen, Scene::CALC_BFIELD);
  scene.render_map (screen);
  screen.write ("test19_bf.bmp");
  maxbf_error = maxef_error = 0;
  scene.calc (screen, Scene::CALC_EFIELD);
  report_error("coil");
  scene.render_map (screen);
  screen.write ("test19_ef.bmp");
  screen.set_axis (vect3d(0, 0, 0), vect3d(1, 1, 0).normal(), vect3d(0, 0, 1));
  scene.calc (screen, Scene::CALC_BFIELD);
  scene.render_map (screen);
  maxbf_error = maxef_error = 0;
  scene.calc (screen, Scene::CALC_EFIELD);
  report_error("coil cross");
  scene.render_map (screen);
  screen.write ("test19_Xef.bmp");
  screen.set_axis (Screen::PLANE_XY);

  do_relative = true;
  scene.calc (screen, Scene::CALC_BFIELD);
  scene.render_map (screen);
  screen.write ("test19_rel_bf.bmp");
  maxbf_error = maxef_error = 0;
  scene.calc (screen, Scene::CALC_EFIELD);
  report_error("rel.coil");
  scene.render_map (screen);
  screen.write ("test19_rel_ef.bmp");
  screen.set_axis (vect3d(0, 0, 0), vect3d(1, 1, 0).normal(), vect3d(0, 0, 1));
  scene.calc (screen, Scene::CALC_BFIELD);
  scene.render_map (screen);
  screen.write ("test19_rel_Xbf.bmp");
  maxbf_error = maxef_error = 0;
  scene.calc (screen, Scene::CALC_EFIELD);
  report_error("rel.coil cross");
  scene.render_map (screen);
  screen.write ("test19_rel_Xef.bmp");
  do_relative = false;
  }

  {
  CoilGen cgen;
  cgen.radius = RADIUS;
  cgen.wr = WR;
  cgen.set_space (SPACING);
  cgen.do_efield = true;
  cgen.shape = CoilGen::ROUND;
  cgen.geom_radial = 4 * RADIUS;
  cgen.geom_axial = 4 * RADIUS;
  cgen.memory_usage = 500;
  cgen.coil_matrix_n = 1;
  cgen.do_cubic = true;
  cgen.generate ("test19_pw.cgen", false, false);

  CoilData cd;
  assert (cd.load ("test19_pw.cgen"));

  Statics s1;
  Statics s2;

  make_polywell_cube (s1, RADIUS, WR, SPACING, 1, 1);
  make_polywell_cube (s2, RADIUS, WR, SPACING, 1, 1, CoilDataFactory<StaticCICU>(cd));

  discard_middle = true;

  Statics statics;
  ElectroDynamics dyn(statics);
  statics.addStaticElement (new DiffStaticElement (s1, s2));
  Scene scene(dyn);
  scene.set_coloring (0, true);

  Screen screen(800, 800, 0.7, 0.7);
  scene.calc (screen, Scene::CALC_BFIELD);
  scene.render_map (screen);
  screen.write ("test19_pw_bf.bmp");
  maxbf_error = maxef_error = 0;
  scene.calc (screen, Scene::CALC_EFIELD);
  report_error("polywell");
  scene.render_map (screen);
  screen.write ("test19_pw_ef.bmp");
  screen.set_axis (vect3d(0, 0, 0), vect3d(1, 1, 0).normal(), vect3d(0, 0, 1));
  scene.calc (screen, Scene::CALC_BFIELD);
  scene.render_map (screen);
  screen.write ("test19_pw_Xbf.bmp");
  maxbf_error = maxef_error = 0;
  scene.calc (screen, Scene::CALC_EFIELD);
  report_error("polywell cross");
  scene.render_map (screen);
  screen.write ("test19_pw_Xef.bmp");
  screen.set_axis (Screen::PLANE_XY);

  do_relative = true;
  scene.calc (screen, Scene::CALC_BFIELD);
  scene.render_map (screen);
  screen.write ("test19_rel_pw_bf.bmp");
  maxbf_error = maxef_error = 0;
  scene.calc (screen, Scene::CALC_EFIELD);
  report_error("rel.polywell");
  scene.render_map (screen);
  screen.write ("test19_rel_pw_ef.bmp");
  screen.set_axis (vect3d(0, 0, 0), vect3d(1, 1, 0).normal(), vect3d(0, 0, 1));
  scene.calc (screen, Scene::CALC_BFIELD);
  scene.render_map (screen);
  screen.write ("test19_rel_pw_Xbf.bmp");
  maxbf_error = maxef_error = 0;
  scene.calc (screen, Scene::CALC_EFIELD);
  report_error("rel.polywell cross");
  scene.render_map (screen);
  screen.write ("test19_rel_pw_Xef.bmp");
  }
}

